﻿using System;
using System.Collections.Generic;
using System.Drawing;

using Emgu.CV;
using Emgu.CV.Structure;
namespace FFTConv
{
    class GaborEmgu
    {
        Matrix<double> KernelRealData;
        Matrix<double> KernelImgData;
        private int GaborWidth = 9, GaborHeight = 9;//滤波器的高度和宽度
        private int orientation = 8, frequency = 3;
        public void Init()
        {
            KernelRealData = new Matrix<double>(GaborWidth, GaborHeight);
            KernelImgData = new Matrix<double>(GaborWidth, GaborHeight);
        }
        #region CalculateKernel && KernelRealPart && KenrelImgPart

        private void CalculateKernel(int Orientation, int Frequency)//计算Gabor核函数
        {
            Init();

            double real, img;

            for (int x = -(GaborWidth - 1) / 2; x < (GaborWidth - 1) / 2 + 1; x++)
                for (int y = -(GaborHeight - 1) / 2; y < (GaborHeight - 1) / 2 + 1; y++)
                {
                    real = KernelRealPart(x, y, Orientation, Frequency);
                    img = KernelImgPart(x, y, Orientation, Frequency);

                    KernelRealData[x + (GaborWidth - 1) / 2, y + (GaborHeight - 1) / 2] = real;
                    KernelImgData[x + (GaborWidth - 1) / 2, y + (GaborHeight - 1) / 2] = img;
                }
        }

        private double KernelRealPart(int x, int y, int Orientation, int Frequency)
        {
            double U, V;
            double Sigma, Kv, Qu;
            double tmp1, tmp2;

            U = Orientation;
            V = Frequency;
            Sigma = 2 * Math.PI * Math.PI;
            Kv = Math.PI * Math.Exp((-(V + 2) / 2) * Math.Log(2, Math.E));
            Qu = U * Math.PI / 8;

            tmp1 = Math.Exp(-(Kv * Kv * (x * x + y * y) / (2 * Sigma)));
            tmp2 = Math.Cos(Kv * Math.Cos(Qu) * x + Kv * Math.Sin(Qu) * y) - Math.Exp(-(Sigma / 2));

            return tmp1 * tmp2 * Kv * Kv / Sigma;
        }

        private double KernelImgPart(int x, int y, int Orientation, int Frequency)
        {
            double U, V;
            double Sigma, Kv, Qu;
            double tmp1, tmp2;

            U = Orientation;
            V = Frequency;
            Sigma = 2 * Math.PI * Math.PI;
            Kv = Math.PI * Math.Exp((-(V + 2) / 2) * Math.Log(2, Math.E));
            Qu = U * Math.PI / 8;

            tmp1 = Math.Exp(-(Kv * Kv * (x * x + y * y) / (2 * Sigma)));
            tmp2 = Math.Sin(Kv * Math.Cos(Qu) * x + Kv * Math.Sin(Qu) * y) - Math.Exp(-(Sigma / 2));

            return tmp1 * tmp2 * Kv * Kv / Sigma;
        }

        #endregion

        public Matrix<double> GaborTransform(Image<Gray, byte> image)
        {
            this.Init();
            for (int m = 0; m < frequency; m++)
                for (int n = 0; n < orientation; n++)
                {
                    this.CalculateKernel(n, m);

                }
                    this.CalculateKernel(0, 0);
            Image<Gray, byte> gabor_real = new Image<Gray, byte>(image.Width, image.Height);
            CvInvoke.cvFilter2D(image.Ptr, gabor_real.Ptr, this.KernelRealData, new Point(-1, -1));
            Image<Gray, byte> gabor_img = new Image<Gray, byte>(image.Width, image.Height);
            CvInvoke.cvFilter2D(image.Ptr, gabor_img.Ptr, this.KernelImgData, new Point(-1, -1));
            Matrix<double> gabor = new Matrix<double>(image.Height, image.Width);
            for (int x = 0; x < image.Width; x++)
                for (int y = 0; y < image.Height; y++)
                {
                    try
                    {
                        double real = gabor_real[y, x].Intensity;
                        double img = gabor_img[y, x].Intensity;
                        gabor[y, x] = Math.Sqrt(real * real + img * img);
                        
                    }
                    catch (Exception exp)
                    {
                        throw exp;
                    }
                }
            return gabor;
            
        }
    }
}
